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Abstract. 

The interstellar medium (ISM) reveals strongly inhomogeneous structures at 
on ■ every scale. These structures do not seem completely random since they obey certain 

on : power laws. Larson's law (1981) a oc R s and the plausible assumption of virial 

equilibrium justify to consider fractals as a possible description of the ISM. In the 
following we investigate how self-gravitation, differential rotation and dissipation 
affect the matter distribution in galaxies. To this end we have performed 3D- 
simulations for self-gravitating local boxes embedded in a larger disk, extending 
the 2D-method of Toomre & Kalnajs (1991) and Wisdom & Tremaine (1988). Our 
\q ■ simulations lead to the conclusion that gravitation, shearing and dissipation can be 

dominantly responsible for maintaining an inhomogeneous and eventually a fractal 
distribution of the matter. 
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1. Introduction 



Keywords: Methods: numerical - Galaxies: ISM - ISM: structure 

ON 
O 

O 

ON 
ON 

The lumpy distribution of matter resulting in Toomre & Kalnajs' (1991, 
hereafter TK) local shearing-sheet experiments, reminds us of the ubiq- 
uitous inhomogeneous state of the ISM as well as the flocculent struc- 
tures of many spirals. Therefore we decided to take up their 2D-model 
and to extend it on 3 dimensions order to investigate more precisely 
the clumpy structure of galactic matter. It is important to include the 
third dimension in the model, because as soon as clumping develops 
dynamical coupling with vertical motion must be strong, contrary to 
the weak coupling existing in a smoothly distributed thin disk. Our 
model considers scales of the order of O(lkpc). Thus we are able to 
investigate the transition regime between the molecular cloud scales 
and the galactic disk scale. Furthermore we can investigate whether 
gravitation and dissipation can maintain the power laws observed in 
molecular clouds on scales at which shearing is dominant. 

This paper is divided up in two main parts. The first part presents 
the model (Sects. 2.1-2.4) and the second part shows the results (Sect. 
3.1) and their analysis, i.e. the determination of the fractal dimension 
(Sect. 3.2) and the verification of Larson's law (Sect. 3.3). 
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2. Model 

2.1. Principle 

In local models only a small part (e.g., everything inside a box with 
a given size) of the system is simulated and more distant regions are 
represented by replicas of the local box. In such a model the orbital mo- 
tion of the particles is determined by Hill's approximation of Newton's 
equations of motion (Hill 1878) 

x - 2Q y = 4fl A x + F x 

y + 2(l x = F y (1) 

z = -u 2 z + F z 

where Aq = —\tq (^r) ro is the Oort constant of differential rotation 
and v is the vertical epicycle frequency. F x ,F y and F z are local forces 
due to the self-gravitating particles. 

Because of the shearing flow the relative positions of the rectangular 
boxes (local box and replicas) change with time, so that a initially 
periodic arrangement of the boxes relative to a fixed Cartesian coor- 
dinate system can not be maintained. As a consequence the forces of 
the self-gravitating particles must be determined by direct summation, 
requiring 0(N 2 ) operations for N particles. 

To increase the performance in our models we calculate the forces 
with the FFT-convolution method, requiring 0(N c \ogN c ) operations, 
where N c is the number of cells, which should have the same order 
of magnitude as the number of particles. This requires a system spa- 
tially periodic at each time. Therefore we use a pair of time-dependent 
affine coordinate systems, whose affinity angle change with time and 
are determined by the shear flow (see Fig. 1). The angle difference 
between the two grids is constant and given by &t&n(L y /L x ). The 
reason to evaluate twice the forces in two different affine grids is for 
avoiding force discontinuities in time when the affine angle would have 
reached a maximum value beyond which a smaller angle would exist. 
After the calculation of the forces for the two affined grids, they are 
transformed into Cartesian coordinates, where they are weighted and 
added. The weighting factors are proportional to the grid inclination 
and normalized. 

2.2. Cooling 

To avoid an increase of the random epicyclic motion due to gravitational 
heating, TK proposed to add artificial cooling forces. Following them 
we include the damping terms — C x x and — C z z in the radial resp. 
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Figure 1. Schematic diagram of the inclined grids seen from above, a.) The initial 
state of the two grids (t = 0). The dashed grid is the forward grid and the solid 
one is the backward grid. Below them, the affinity angles of the grids are indicated, 
ctf resp. ah- The angle between the two grids remains the same for all times, b. ) 
The grids at t = L y /(2L x Qo), where a/ = — at- c.) When the grids attain these 
inclinations, they jump back to the positions shown in (a) and the process starts 
again without discontinuity in the dynamics. 



vertical forces (F x , F z ) controlling the particle motions via Eq. (1). Two 
different damping terms are necessary to reflect the different directional 
collision rates induced by the anisotropic velocity ellipsoid. The cooling 
coefficients are chosen in such a way that the velocity dispersions of the 
random motions reach a stable level and don't differ much from their 
initial values. The same stability conditions holds for the disk scale 
height zq. 
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Figure 2. The evolution of the particle positions seen from above. The number of 
rotations of the shearing box around the galactic center are indicated above each 
panel. The number density at the beginning of the simulation is n = lOlOOA^, 
corresponding to 32724 particles. 



3. Results 

3.1. Shearing boxes 

Fig. 2* and 3 reveals the spatial distribution of matter in a box, rep- 
resenting a section of a galactic disk. The size of the box sides [L x x 
L y x L z ) is given by the critical wavelength A cr it , defining the scale for 
which the theory of swing amplification predicts the strongest response 
(Toomre 1981). 

Because we are interested in long time behavior, the simulations 
were performed for t = 20 galactic rotations. In the initial state t = 
the particles are distributed uniformly in the x-y-plane. The particle 
distribution in z-direction obeys p oc sech 2 (2;/zo), where p is the density 
and zq is the disk scale height. The velocities at t = are determined 
by the shear-flow 

x = 0, y = -2A x, z = (2) 

and the Schwarzschild velocity ellipsoid. 

* Full resolution paper available at http://obswww.unige.ch/Preprints/cgi- 
bin/Preprintshtml.cgi?#DYNAMIC. 
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Figure 3. The positions of the particles inside the slice, —0.1 A cr it < y < 0.1 A cr it, 
seen in the direction of orbital motion. 



Like the 2D-model of TK, our 3D-model reveals a fast "structure 
formation" . After the first rotations the striations are already developed 
and are more or less maintained during the rest of the simulation. 

In our model t osc < t C ooi,a; < t coo \ iZ is valid, where t osc is the period 
of the unforced epicyclic motion. t coo \ jX and t coo \ jZ are the cooling times 
for the radial resp. vertical damping in fact t coo i,x °c C~ l and i coo /,z oc 
Cy 1 . Furthermore, the cooling times depend on the particle number 
density n. For the simulation in Fig. 2 the following relation is valid: 

^osc • icoo\,x • ^cool,z ~ 1 • 40 . 300. 

The random velocity dispersions and the disk scale height were cal- 
culated after every half rotation and are plotted in Fig. 4. Gravitational 
instability and shearing, which are responsible for an important part 
of the striations, have an effect particularly in the x-y-plane and leads 
to a fast increase of a x and a y . 

3.2. Fractal dimension 

During this studies our interest was particularly focused on the follow- 
ing question: Can self-gravitation, shearing and dissipation be domi- 
nantly responsible for maintaining a fractal distribution of the matter 
in a disk and can they account for Larson's law (Larson 1981)? To 
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Figure 4- The left panel shows the velocity dispersions a x , a y and a z as a function 
of time t, indicated in galactic rotation unit. The right panel reveals the evolution 
of the disk scale height zq. 



answer this question we calculated the appropriate scaling laws and 
determined the fractal dimension for the structures shown in Fig. 2 
and 3. 

We determined the fractal dimension, quantifying how much space 
our system fills, as follows: We choose a representative set of particles 
and count for each particle the number of neighboring particles N(R) 
inside a certain radius R. If we repeat this for other values of R we can 
find the fractal dimension D via 

N(R) oc R D . (3) 

Because the structures we examine are not an idealized mathematical 
set, but model a physical system, we have to take into account upper 
and lower cutoffs in the analysis. An upper limit due to the numerical 
model is given by the size of the simulation box in the x-y-plane. On this 
scale the system becomes periodic, meaning that it can not be fractal 
on scales 0(L X ). A lower limit is due to the finite resolution of the box 
grid. If the grid cells have the size l x x l y x l z and I = l x = l y > l z then 
we can't expect fractal structures below 21. Fig. 5 reveals the fractal 
dimension D = dln(N) / 'dln(R) as a function of the radius R, i.e. of the 
scale. The solid line corresponds to the initial distribution of matter: 
p(x,y) = const., p(z) oc sech 2 (z/zo)- The other lines give the mean 
dimensions of the structures for the simulation terminal phase. The 
intensity of the structures and thus the structure dimension D depend 
on the cooling. To show this clearly we performed simulations with 
varying cooling forces. The figure reveals clearly that the stronger the 
cooling, the more intensive the structures and the lower the structure 
dimension D are. The structures in Fig. 2 and 3 were performed with 
a "strong" cooling. For this structures D runs in a thin band between 
0.2A cr i t and the upper model limit with a minimum of D ps 1.83 at 
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Figure 5. The structure dimension D as a function of the scale R. The solid 
line corresponds to the initial state. On large scales this state represents a 2D 
matter distribution, whereas on small scales 7? < zo it tends to D > 2.6. The 
other curves represent the simulation terminal phase for different cooling forces 
(weak, middle, strong). The relative magnitudes of the cooling coefficients are: 
Ca:,weak '• Cx, middle '• Cx, strong ?s 1 : 1.5 : 2. Above the figure the natural logarithm of 
the lower model limit is indicated. 



R = 0.35A cr it- Pfenniger (1996) showed, that a fractal-like, hierarchical 
gravitating system in statistical equilibrium can only exist with D < 2. 
In the range, in which D runs in a thin band this condition is fulfilled, 
but the dynamical range is too small to call the structures fractals. 



3.3. Larson's law 



If the fractal organization of the matter extends to scales at which the 
shear flow becomes important and the system is still virialized, then 
the validity range of Larson's law 

<t cc R or 5 = — (4) 

dlnR w 

should exceed the size of Giant Molecular Clouds. We checked Larson's 
law for scales > 100 pc by analyzing the phase-space of the particles in 
Fig. 2. 

Fig. 6 shows 5 as a function of the scale R. There is no plateau 
between 0.3 < 5 < 0.5, but 8 reaches this range for the bigger scales, 
where also 25 + 1 = D is roughly fulfilled. 
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Figure 6. 8 of Larson's law as a function of the scale R. To calculate 8 we used the 
simulation with the particle density n = 40450 and the resolution I = 0.03A crit The 
solid line corresponds to the initial state, whereas the dashed line indicate 8 for the 
simulation terminal phase. 



4. Conclusion 

Our local 3D-simulations show, that self-gravitation and dissipation 
ensure a statistical equilibrium at scales at which the shear flow is 
important. Moreover they reveal a fractal-like distribution of matter 
on kpc-scales. The fractal dimension found in the simulations can be 
lower than 2, as observed in molecular clouds. The specific value found 
depends on the relative strengths of the competing gravitational and 
dissipation processes. Because the dynamical range of our simulations 
is still small, it would be premature to call the found structures strict 
fractals. The anisotropy of the velocity-dispersion ellipsoid, resulting 
from our simulations, has systematically the same ordering (ctr > ct$ > 
a z ) as observed in the Galaxy and in N-body simulations of spirals. 
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